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Superconductivity in the Two-Dimensional t-J Model at Low 

Hole Doping 

E. S. Heeb and T. M. Rice 
Abstract 

By combining a generalized Lanczos scheme with the variational Monte Carlo 
method we can optimize the short- and long-range properties of the ground- 
state separately. This allows us to measure the long-range order of the ground- 
state of the t-J model as a function of the coupling constant J/t, and identify 
a region of finite d-wave superconducting long-range order. With a lattice 

size of 50 sites we can reliably examine hole densities down to 0.16. 
PACS numbers: 71.10,+x, 71.27.+a, 74.20.Mn 
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Presently one of the most interesting questions in the study of strongly correlated elec- 
tron systems is to determine the region of superconductivity in the phase diagram of the 
two-dimensional t-J model Despite considerable effort, high temperature expansion 

0fj] and Quantum Monte Carlo calculations |5],|6[] are unable to provide conclusive evidence 
for or against superconductivity. Although variational calculations ]7| are in favor of super- 
conductivity, they are not able to establish independent evidence in a region of the phase 
diagram where magnetically ordered phases are competing with superconductivity and have 
almost the same energy. To date the strongest unbiased indication for superconducting order 
in the groundstate comes from exact diagonalization of a 4 x 4 cluster || . With such a small 
cluster only a few filling fractions are available and it remains open how much the finite size 
effects contribute to the results. In this letter we will present results on VoO x v^O cluster, 
which allows us to reach smaller hole dopings and reduce the finite size effects. By varying 
the long-range and short-range correlations separately we can go beyond a standard vari- 
ational approach and we are able to measure the long-range pair-pair correlation function. 
This enables us to identify a region of a d-wave superconducting groundstate in the phase 
diagram. 

The t-J model has attracted considerable attention, because it is the simplest model to 
describe strongly correlated electrons. The perturbative approaches are not possible since 
there is no exact solution that could be used as a starting point. In contrast, many of 
the numerical methods that have been used for strongly correlated electrons do not rely on 
a small parameter in the Hamiltonian. Exact diagonalization provides reliable results of a 
variety of properties. However, for the two-dimensional t-J model only systems up to a size of 
4x4 have been investigated for all fillings . The results for bigger systems are restricted to 



certain fillings only, e.g., 2 holes on 26 sites fllQfl. The limited number of systems and fillings 



makes it difficult, if not impossible, to obtain reliable information about the thermodynamic 
limit. Properties that are related to the short range behavior like the energy (e.g., total 
energy, binding energy, spectral function) are more reliable than those involving long-range 
behavior. Indeed, the evidence for the presence or absence of superconducting long-range 



order from exact diagonalization is very limited |||ll[] . High-temperature expansions directly 
lead to results which are valid in the thermodynamic limit but the extrapolation from 
finite temperatures to T = is difficult. So far, only the equal time behavior of the spin 
and charge degrees of freedom have been successfully analyzed while the question of 
superconducting order remains open. The Quantum Monte Carlo methods |T2[ which are 



very powerful for weakly interacting systems show severe restrictions due to the fermion 
sign problem especially in the strong coupling limit. The available results cover mainly the 
intermediate coupling regime of the Hubbard model f| whereas for the inherently strong 
coupling t-J model only the cases of one and two holes have been considered ||. 

With variational approaches no fermion sign problem occurs and the systems are big 
enough to show only small finite size effects. If enough information about the symmetry 
of the groundstate is known, a variational wavefunction can be constructed to model this 
groundstate. Typically, such a variational wavefunction is given analytically from a mean 
field ansatz, so that even for fermions expectation values can be readily evaluated by Monte 
Carlo sampling ||13||. It is well known that in the t-J model various phases with different 



broken symmetries compete with each other. The energies of these phases will be close 
to each other at small dopings so that the energy differences are comparable to the error 
introduced by using modified mean field forms for the variational wavef unctions. The results 
of these variational studies remain therefore inconclusive to some extent, and the regions of 
stability of the various phases in the phase diagram can only be estimated qualitatively 0. 
Systematic iterative improvements have been used to remove the bias in the choice of 



the wavefunction. These methods range from the power method to Lanczos iterations 
H| . While they do remove the bias, these methods are restricted to a few iterations only. 
The computing time needed to reduce the statistical error increases rapidly with the number 
of iterations. In this letter we will use a generalized Lanczos approach which optimizes the 
short- and long-range correlations separately and extract the relevant information from one 
iteration. 

The Hamiltonian for the t-J model is defined in the subspace with no doubly occupied 



sites as 

H = ~t ( 6 L 5 J> + h.c.) + J ( S i- S i - \ n i n j) (!) 

<i,j>,(T <i,3> 

where c\ a = c\ a (1 — ni- a ) prevent doubly occupied sites and the rest of the notation 
is standard. We perform our calculations on a 50 sites cluster with periodic boundary 
conditions and periods (7, 1) and (—1,7). Exact diagonalization || has found that at the 
largest distance available in the 4x4 cluster the pair-pair correlation function C(R) = 
(l/N)J2i< A|A i+ R > corresponding to a nearest-neighbor singlet d-wave pairing operator 
Aj = (1/2) EaQ-a {ci+x,a + c%-x,a - Ci+y,a - (k-s,o) is dominant. This long-range order is 
realized in a mean field wavefunction as follows 

|* (D, rf) = V G V N J] («k + vAA*,l) 1°) ( 2 ) 

k 

where Vq and Pn are the Gutzwiller- and N-particle projectors, respectively. The ratio 
^k/^k = A. k / (l^k + y^l~+ with ^ k = —2 (cos (k x ) + cos (k y )) — jjl has the standard BCS 
form. With A k = D (cos (k^) — cos(k 3/ )) this wavefunction has by construction a finite 
long-range d-wave order in the thermodynamic limit for a finite D-parameter, whereas for 
D = it reduces to the Gutzwiller projected Fermi sea. For the k points where A k has a 
node the ratio v^/u^ is not well defined when £ k < 0. In the thermodynamic limit the nodes 
of A k in the Brillouin zone are negligible. The effect of these nodes on the wavefunction 
accounts for much of the finite size effects. Due to the tilted periodic boundary conditions 
the 50 sites lattice has only one point with A k = (at k = 0) and is thus an optimal choice 
to reduce the finite size effects. Since k = is deep inside the Fermi sea we set v k =o — » 1 and 
« k=0 — > which leads to a ratio f k /w k — > oo. In an actual calculation we choose a large but 
finite ratio. For D — » the choice of this ratio has a bigger influence on the wavefunction 
and the Fermi sea will be defined as the extrapolation from small but finite values of D. It 
is important to note that this wavefunction is constructed to display a specific long-range 
behavior and that there is no direct control over the short-range part. The Hamiltonian with 
its nearest-neighbor terms may therefore well favor a gap parameter D which is shifted away 
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from the value that would correspond to the correct long-range behavior of the groundstate. 
Standard variational calculations provide no control that would allow one to find out whether 
D, which determines the long-range correlations, is over- or underestimated. This is one of 
the main disadvantages of the standard variational calculations. 

We remedy this situation by optimizing the short-range correlations independently from 
the variational wavefunction. Since the Hamiltonian only contains nearest-neighbor terms, 
we construct the most general nearest-neighbor operator, which conserves the quantum 
numbers for the spin and space symmetries: 



In the combined wavefunction A (D)) we can adjust for the best short-range correlations 
through the choice of the parameters The parameter D now only controls the long-range 
behavior for which it was designed. By using an operator A with the same length scale as 
the Hamiltonian 7i we arrive at a scheme, which is similar to a Lanczos iteration. However, 
we allow the parameters ctj to be adjusted independently of the coupling strengths in the 
Hamiltonian. Our approach can therefore be regarded as a generalization of the Lanczos 
scheme. 

In standard variational calculations the Rayleigh-Ritz principle (RR) is used to find 
the best variational parameters. The expectation value E RR = min^, H / (^l^) is 
the lowest upper bound for the groundstate energy that can be achieved with a given set 
of variational wavefunctions \^/). Additionally the variance = (H 2 ) — (H) 2 is small 



if a wavefunction is close to an eigenstate of the Hamiltonian. When approaching the 
groundstate, both of these values should become smaller. Although, this is a necessary 
condition, it is not sufficient to distinguish from the case where one approaches an excited 
eigenstate. 

In our approach we use the same Rayleigh-Ritz principle but with the addition of the 
generalized Lanczos operators (RRGL). We arrive at the expectation value 



A = a + Oil J2 ( C t,a C j,v + h.c.) + "2 J2 S i " S i + a 3 



(3) 
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which is still a lowest upper bound for the groundstate energy. Also the variance is again 
used to estimate the width of the energy spectrum of the wavefunction. Additionally to the 
standard variational approach we can now compare the RR- to the RRGL-wavefunction. 
This allows us to judge whether the mean field parameter D in RR is indeed shifted to adjust 
for the short-range behavior or whether the long-range correlations are maintained in the 
groundstate. Furthermore, we know from the Lanczos method used in exact diagonalization 
that when one starts from a state with an energy spectrum which is centered around an 
excited eigenstate, the first iteration will already redistribute the weights towards lower 
energies, such that the variance will first increase. In this way we use the variance as a 
unbiased indicator to judge the quality of the RR- wavefunction. Using these criteria the 
wavefunctions Eq. fl2|), with d-wave order parameter, prove to be a good choice consistent 
with the results from exact diagonalization ||. In this work we will therefore concentrate 
on the set of wavefunctions Eq. (|2]) as our starting point. 

In Fig. [1] we show the energy for the RR wavefunction and the improved RRGL wave- 
function at quarter filling. For 16 sites this corresponds to 8 holes whereas for 50 sites we use 
24 holes, so as to have an even number of particles. For the case of 8 holes on 16 sites we can 
compare our results to the exact groundstate energy and we find that the short-range corre- 
lations account for about 80 % of the missing correlation energy. Fig. ^| shows the square root 
of the variance, i.e. the width of the energy spectrum. This value is considerably reduced for 
the RRGL wavefunction consistent with a wavefunction which is closer to the groundstate. 
For the 50 sites lattice there are no results for the exact groundstate energy available and we 
have to estimate how the various quantities scale with system size. While the wavefunction 
scales with the size of the system, the generalized Lanczos operators A always act on the 
same length scale as the Hamiltonian. The energy and will therefore scale with the 
system size, a-n is also the energy scale for the improvement A E = E RK — -Errgl- If the 
wavefunction has the right long-range behavior then the operators A need to improve the 
correlations only on the same length scale as the Hamiltonian and the ratio Ae/ct-h should 
remain constant. Indeed we find this value to be slightly bigger for the 50 sites lattice than 
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for 16 sites. This again supports the observation that the wavefunction Eq. (Q) describes 
the correct long-range behavior of the groundstate. For the 50 sites lattice we investigated 
several fillings, which all showed analogous results to the ones described above. Specifically 
we looked at the closed shell configurations of 8, 16, and 24 holes. 

To investigate the superconductivity we measure the pair-pair correlation function C(R). 
We find that for the 50 sites lattice C(R) is flat for the larger distances indicating that the 
finite size effects are small. We can therefore take Coo = C(-R max ) as a measure for long-range 
order. In the standard variational approach Coo is a monotonic function of D and contains 
no additional information. With our new method we can now test how the introduction of 
the operators A in the wavefunction affects C^. If we start with too much long-range order 
the operators A will redistribute the weight in the correlation function and suppress C^. On 
the other hand too small a value for Coo will be enhanced. If we start with the correct long- 
range order that corresponds to the groundstate, the operators A will only affect the short 
range part of C(R) and will be unaffected. In that case we have effectively separated 
the short- and long-range parts of the wavefunction. 

We will illustrate this for the case of 8 holes on 50 sites corresponding to a hole density 
of 0.16. In Fig. 0(a) we show C(R) for one value of D = OAt and /i = — 0.8t. The solid line 
corresponds to the RR-wavefunction. We can see that the long-range tail is well saturated. 
For J < J c the long-range correlations are suppressed, while for J > J c they are enhanced. 
This is shown by the dashed lines. For D = OAt we find J c ~ l.Ot. We can now combine the 
data obtained for different gap parameters D. In Fig. |3](b) we show Coo as a function of the 
coupling constant J/t. The solid line again corresponds to the RR- value while the dashed 
line shows the suppression and enhancement for the RRGL- values. The point where the solid 
and dashed lines cross is the value Coo which remains unchanged under iteration and we 
take this as the long-range order Coo of the groundstate. For other values of the variational 
parameter D we repeat the same procedure and obtain the points shown in Fig. 0(b). The 
error bars in J/t indicate the region where the suppression or enhancement is within one 
standard deviation. We can thus map out C^J/t) for the groundstate and determine the 
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critical value J c = (0.44 ± 0.04)t for the onset of superconducting long-range order. For the 
other closed shell configurations of 16 and 24 holes we obtain analogously the critical values 
J c = (0.62 ± 0.08)t and J c = (0.78 ± 0.11)t, respectively. Using these points we construct 
the phase diagram Fig. |j. Ohta et al. JTIJ] observe a finite superconducting gap also at lower 
values of J/t for 16 and 18 sites systems. However, they only report results on open shells, 
in order to have k points at the Fermi surface. For open shells the Fermi sea is degenerate 
and the d-wave condensate might be favored in order to remove this degeneracy. Such finite 
size effects are reduced for the closed shell configurations used in this work. 

The wavefunctions used in this work describe a homogeneous electron distribution for 
all variational parameters. For large values of J, the t-J model exhibits phase separation 
T6|,|3[|. In that region of the phase diagram, the groundstate is not well represented by the 



homogeneous wavefunctions. However, since the phase separated state is a mixture of two 
homogeneous states — the half filled Heisenberg antiferromagnet (hole doping 5 = 0) and 
a state with finite hole doping — we can use the Maxwell construction to obtain its energy. 
For any fixed value of J the energy of the homogeneous wavefunctions will be a smooth 
function of the hole density 5, which we can describe by a polynomial in 5. In the region of 
phase separation this polynomial will be curved downwards, so that a combination of two 
homogeneous states (represented by a straight line) will lower the energy. The maximal hole 
density S C (J), which lowers the energy between and 5 C determines the phase separation 
line as a function of J. This line is also shown in Fig. |j. For the Maxwell construction we 
only use the energies of the closed shell configurations (8, 16, and 24 holes), for which we 
can expect the finite size effects to be minimal. We can then identify the region of d-wave 
superconducting long-range order from close to quarter filling (5 = 0.48) down to a hole 
density of 5 = 0.16. This is shown by the shaded region in the phase diagram Fig. [|. 

In conclusion we have presented a method which allows us to measure the long-range 
behavior of the groundstate by separating long-range from the short-range contributions. 
This allows us to calculate the long-range d-wave pair-pair correlation function, which is 
a direct measure for the superconducting order parameter. We identify the region of d- 
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wave superconducting long-range order down to a hole density of 5 = 0.16. For smaller 
hole densities the calculations will have to be extended to larger lattice sizes. Our method 
extends the results from the exact diagonalization, as it can be applied to larger systems. 

We would like to acknowledge useful discussions with S. Barnes, E. Dagotto, G. Blatter, 
N. Bonesteel, C. Gros, R. Joynt, B. Koltenbah, T. K. Lee, H. Tsunetsugu, and F. C. Zhang. 
We thank the Swiss Nationalfond for financial support. 
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FIGURES 

FIG. 1. Energy per site in units of t as a function of coupling constant J/t for the t-J model 
at quarter filling. .Err is the energy for the standard variational (Rayleigh-Ritz, RR) approach, 
whereas Errgl shows the energy after optimizing the short-range correlations with generalized 
Lanczos operators. For 16 sites the values of the exact diagonalization from Dagotto et al. (dotted 
line) are shown for comparison. The lines are upper bounds for the groundstate energy. The 
statistical error is smaller than the line width. 

FIG. 2. Variance as a function of coupling constant J/t. Shown is the square root £pft of the 
variance per site in units of t. This value measures the width of the energy spectrum and is a 
measure for the quality of the wavefunction. RR and RRGL correspond to the same wavefunctions 
as in Fig. 1 

FIG. 3. (a) d-wave pair-pair correlation function C(R) as a function of distance for a hole den- 
sity of 0.16 and a gap parameter D = OAt. The solid line corresponds to the raw RR wavefunction. 
The dashed lines show the long range correlation for the RRGL improvement for a value of J > Jo 
and J < Jq respectively. At the critical value Jo the long-range correlation is unchanged from the 
RR ansatz. (b) long-range d-wave correlation as a function of coupling constant J/t. The solid 
and dashed lines show for the same gap parameter as in (a). For other variational parameters 
only the points Coo(J/t) where is unchanged are shown. 

FIG. 4. Phase diagram of the t-J model with parameters hole density 5 and coupling con- 
stant J/t. The solid line shows the line of phase separation. The points indicate the onset of 
super conduct vity. In the shaded region we observe a finite superconducting long-range order. 
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